

cuts<-c("all","married","prerich","ownhome")
cutnames<-c("_agg","","rich","_homeowner")
cutlabs<-list(c(""),
              c("Single","Partnered"),
              c("Low-Wealth","High-Wealth"),
              c("Renter","Owner"))

cutlims<-list(c(0,0.5),
              c(0,0.5),
              c(0,0.6),
              c(0,0.5))

for(cutnum in 1:length(cuts)){
  cut<-cuts[cutnum]
  cutname<-cutnames[cutnum]
  cutlab<-cutlabs[[cutnum]]
  cutlim<-cutlims[[cutnum]]

  
    linetype<-c("dashed","solid")
    vref<-1

  eventmeans<-TRUE
  source("./psid_eventpreamble.R")
  
  
    rawdata<-rawdata[educ<=2,]
    rawdata_mod<-rawdata_mod[educ<=2,]

  #Stats on worker's comp vs. SS:
  rawdata_mod[,everhealth:=max(health),by=newid]
  rawdata_mod[,evercomp:=max(wcomp>0,na.rm=TRUE),by=newid]
  rawdata_mod[,everss:=max(soc_sec>0,na.rm=TRUE),by=newid]
  rawdata_mod[,.(mean(soc_sec),mean(wcomp,na.rm=TRUE),.N),by=.(health,marnow)]
  rawdata_mod[age==onset_age,.(mean(everss),mean(evercomp),.N),by=.(everhealth)]

meandata<-rawdata[time>= -6& time  <= 20 & empyet == 1,.(anyDI=mean(anyDI),N=sum(N)),by=c("time",cut)]
meandata<-meandata[order(time),]


meandata_mod<-rawdata_mod[time>= -6& time  <= 20 & empyet == 1 &!is.na(anyDI),.(anyDI=mean(anyDI),N=sum(N)),by=c("time",cut)]
meandata_mod<-meandata_mod[order(time),]

meandata<-meandata[meandata_mod,,on=c("time",cut)]
setnames(meandata,c("i.anyDI","i.N"),c("anyDI_mod","N_mod"))
pdf(paste0("./PSID_eventstudy/headsev_level_anyDI",cutname,".pdf"))
print(ggplot(data= meandata[time>=-4 & time <=4 & !is.na(get(cut)),])+
  geom_line(aes(x=as.numeric(as.character(time)),y=anyDI,group=as.factor(get(cut)),linetype = as.factor(get(cut))),size=1)
+geom_vline(xintercept=0,linetype="dashed")
+geom_hline(yintercept=0,linetype="dashed")
+ylim(cutlim)
+scale_x_continuous(breaks=-4:4)
+scale_linetype_manual(labels=cutlab,values=linetype)
+theme(legend.position="bottom")
+labs(x="Years from Onset",y="Share Receiving DI Benefits",linetype="")
+theme(legend.key.size = unit(2,"line"))
)
dev.off()

rawmeans<-feols(as.formula(paste0("anyDI ~ as.factor(",cut,")-1")), data = rawdata[time>= 0 & time <= 4 & empyet==1,],cluster="newid")$coeftable[,c("Estimate","Std. Error")]
write.csv(rawmeans,file=paste0("./PSID_eventstudy/",paste(c("sev",cutname,"_eventmeans_table.csv"),collapse="")), 
          row.names = TRUE, col.names = FALSE)

if(cut!="all"){
  vvals <- unique(rawdata[,get(cut)])
  vvals<-vvals[!is.na(vvals) & vvals!=vref]
  
rawdata[,agebin:=floor(age/10)*10]
rawdata[,agebin:=max(agebin * (time==0)),by=newid]

rawdata[,pweight:=NULL]
rawdata[,dummy:=get(cut)==vref]
for(vlev in vvals){
  rawdata[,eval(paste0("pval_",vlev)):=NULL]
  rawdata[!is.na(agebin) & get(cut)%in%c(vref,vlev) 
            ,eval(paste0("pval_",vlev)):=feols(as.formula(paste0("dummy~ 1 | interaction(time,agebin)")),
                                               data = rawdata[!is.na(agebin)  & get(cut)%in%c(vref,vlev),
                                               ],)$fitted.values]
  rawdata[get(cut)==vlev & get(paste0("pval_",vlev)) < 1 & get(paste0("pval_",vlev)) > 0,
            pweight:=get(paste0("pval_",vlev))/(1-get(paste0("pval_",vlev)))]
}
#Imposing common support:
for(vlev in vvals){
  rawdata[get(cut)==vref & !is.na(get(paste0("pval_",vlev))),pweight:=1]
  rawdata[(round(get(paste0("pval_",vlev)),5) == 1 | round(get(paste0("pval_",vlev)),5) == 0),pweight:=NA]
}


rawdata[time==0&!is.na(pweight),weighted.mean(age,pweight),by=c(cut)]
rawdata[time==0,mean(age),by=c(cut)]

write.csv(feols(as.formula(paste0("anyDI ~ as.factor(time)+as.factor(time):",cut,"-1")), data = rawdata[time>= 0 & time <= 4 & empyet==1,],cluster="newid",
      weights=rawdata[time>= 0 & time <= 4 & empyet==1,]$pweight)$coeftable,
      file=paste0("./PSID_eventstudy/",paste(c("htest_sev_dynamic",cutname,".csv"),collapse="")), 
      row.names = TRUE, col.names = FALSE)

write.csv(feols(as.formula(paste0("anyDI ~ as.factor(",cut,")")), data = rawdata[time>= 0 & time <= 4 & empyet==1,],cluster="newid",
                weights=rawdata[time>= 0 & time <= 4 & empyet==1,]$pweight)$coeftable,
          file=paste0("./PSID_eventstudy/",paste(c("htest_sev_pooled",cutname,".csv"),collapse="")), 
          row.names = TRUE, col.names = FALSE)


meandata<-rawdata[time>= -6& time  <= 20 & empyet == 1 &!is.na(pweight),.(anyDI=weighted.mean(anyDI,pweight),N=sum(N)),by=c("time",cut)]
meandata<-meandata[order(time),]


pdf(paste0("./PSID_eventstudy/headsev_agebal_level_anyDI",cutname,".pdf"))
print(ggplot(data= meandata[time>=-4 & time <=4,])+
        geom_line(aes(x=as.numeric(as.character(time)),y=anyDI,group=as.factor(get(cut)),linetype = as.factor(get(cut))),size=1)
      +geom_vline(xintercept=0,linetype="dashed")
      +geom_hline(yintercept=0,linetype="dashed")
      +ylim(cutlim)
      +scale_x_continuous(breaks=-4:4)
      +scale_linetype_manual(labels=cutlab,values=linetype)
      +theme(legend.position="bottom")
      +labs(x="Years from Onset",y="Share Receiving DI Benefits",linetype="")
      +theme(legend.key.size = unit(2,"line"))
)
dev.off()
}

pdf(paste0("./PSID_eventstudy/headmod_level_anyDI",cutname,".pdf"))
print(ggplot(data= meandata_mod[time>=-4 & time <=4,])+
        geom_line(aes(x=as.numeric(as.character(time)),y=anyDI,group=as.factor(get(cut)),linetype = as.factor(get(cut))),size=1)
      +geom_vline(xintercept=0,linetype="dashed")
      +geom_hline(yintercept=0,linetype="dashed")
      +ylim(cutlim)
      +scale_x_continuous(breaks=-4:4)
      +scale_linetype_manual(labels=cutlab,values=c("dashed","solid"))
      +theme(legend.position="bottom")
      +labs(x="Years from Onset",y="Share Receiving DI Benefits",linetype="")
      +theme(legend.key.size = unit(2,"line"))
)
dev.off()

rawmeans<-feols(as.formula(paste0("anyDI ~ as.factor(",cut,")-1")), data = rawdata_mod[time>= 0 & time <= 4 & empyet==1,],cluster="newid")$coeftable[,c("Estimate","Std. Error")]
write.csv(rawmeans,file=paste0("./PSID_eventstudy/",paste(c("mod",cutname,"_eventmeans_table.csv"),collapse="")), 
          row.names = TRUE, col.names = FALSE)

if(cut != "all"){
rawdata_mod[,agebin:=floor(age/10)*10]
rawdata_mod[,agebin:=max(agebin * (time==0)),by=newid]
rawdata_mod[!is.na(get(cut)),pval:=feols(as.formula(paste0(cut," ~ 1 | interaction(time,agebin)")),
                     data = rawdata_mod)$fitted.values]

rawdata_mod[,pweight:=NULL]
rawdata_mod[get(cut)==1 & pval < 1 & pval > 0,pweight:=1]
rawdata_mod[get(cut)==0 & pval < 1 & pval > 0,pweight:=pval/(1-pval)]
rawdata_mod[time==0&!is.na(pweight),weighted.mean(age,pweight),by=c(cut)]
rawdata_mod[time==0,mean(age),by=c(cut)]



meandata_mod<-rawdata_mod[time>= -6& time  <= 20 & empyet == 1 &!is.na(pweight) &!is.na(anyDI),.(anyDI=weighted.mean(anyDI,pweight),N=sum(N)),by=c("time",cut)]
meandata_mod<-meandata_mod[order(time),]



write.csv(feols(as.formula(paste0("anyDI ~ as.factor(time)+as.factor(time):",cut,"-1")), data = rawdata_mod[time>= 0 & time <= 4 & empyet==1,],cluster="newid",
                weights=rawdata_mod[time>= 0 & time <= 4 & empyet==1,]$pweight)$coeftable,
          file=paste0("./PSID_eventstudy/",paste(c("htest_mod_dynamic",cutname,".csv"),collapse="")), 
          row.names = TRUE, col.names = FALSE)

write.csv(feols(as.formula(paste0("anyDI ~ ",cut)), data = rawdata_mod[time>= 0 & time <= 4 & empyet==1,],cluster="newid",
                weights=rawdata_mod[time>= 0 & time <= 4 & empyet==1,]$pweight)$coeftable,
          file=paste0("./PSID_eventstudy/",paste(c("htest_mod_pooled",cutname,".csv"),collapse="")), 
          row.names = TRUE, col.names = FALSE)

pdf(paste0("./PSID_eventstudy/headmod_agebal_level_anyDI",cutname,".pdf"))
print(ggplot(data= meandata_mod[time>=-4 & time <=4,])+
        geom_line(aes(x=as.numeric(as.character(time)),y=anyDI,group=as.factor(get(cut)),linetype = as.factor(get(cut))),size=1)
      +geom_vline(xintercept=0,linetype="dashed")
      +geom_hline(yintercept=0,linetype="dashed")
      +ylim(cutlim)
      +scale_x_continuous(breaks=-4:4)
      +scale_linetype_manual(labels=cutlab,values=c("dashed","solid"))
      +theme(legend.position="bottom")
      +labs(x="Years from Onset",y="Share Receiving DI Benefits",linetype="")
      +theme(legend.key.size = unit(2,"line"))
)
dev.off()
}
}
